#include<stdio.h>
#include<stdlib.h>
#include <mpi.h>
#include <omp.h>
#include <string.h>

float	**Alocar_matriz_real (int, int);
float	**Liberar_matriz_real (int, int, float **);
float	*Alocar_vetor_real (int);
float	*Liberar_vetor_real (int, float *);
void	Le_matriz_vetor(int,float *, float **,FILE *);
float	Numero_maximo(int,float *);
int	*Alocar_vetor_int(int );

#define BLOCK_LOW(id,p,n) (((id)*(n))/(p))
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,(p),(n)) - 1)
#define NTHREAD 2

int main(int argc, char *argv[])
{

	int	ordemMatriz;			// J_ORDER
	int	filaAvaliacao;			// J_ROW_TEST
	float	erro,mr,valor_filaAvaliacao;	// J_ERROR
	int	interacoes;			// J_ITE_MAX
	float	**matriz;
	float	*vetor,*vetor_filaAvaliacao;
	float	*result_pos,*result_antes,*dif;
	int	i,j,k,l;				//Variaveis de manipulacao das matrizes
	float	aux;				//variavel auxiliar
	int	nl,posicao;			// numero de linhas e posicao reponsavel na matiz
	int	flag,resto;
	int 	*rcv_disp, *rcv_cnt;		//variaveis para MPI_Allgetherv
	float	*x,*xg;
	FILE *fp;


	/*Variáveis MPI */   
        int             dest,noProcesses, processId, src, tag;
        MPI_Status      status;   

        /* Inicializacção do MPI */   
        MPI_Init(&argc, &argv);   
        MPI_Comm_size(MPI_COMM_WORLD, &noProcesses);
        MPI_Comm_rank(MPI_COMM_WORLD, &processId);

	//abre o arquivo para leitura
	if ( (fp = fopen(argv[1],"r")) == NULL ) 
        {
                printf("Erro ao abrir o arquivo.\n");
                exit(1);
        }

	//leitura dos dados iniciais
	fscanf(fp,"%d",&ordemMatriz);
	fscanf(fp,"%d",&filaAvaliacao);
	fscanf(fp,"%f",&erro);
	fscanf(fp,"%d",&interacoes);

	 //divide o numero de linha que cada processo deve ler
	resto = ordemMatriz%noProcesses;
        nl=(int)(ordemMatriz/noProcesses);
        if(processId==0) 
	{
		nl=nl+resto;
		posicao=0;
	}
	else
	{
		posicao=(nl*processId)+resto;
	}
	
	//aloca na memoria a matriz e o Vetor
	xg = Alocar_vetor_real(ordemMatriz);
	x = Alocar_vetor_real(nl);// resultado parcial de cada nó
	matriz = Alocar_matriz_real (ordemMatriz, ordemMatriz); // matriz quadrada
	vetor = Alocar_vetor_real (ordemMatriz);
	result_pos = Alocar_vetor_real (ordemMatriz);
	result_antes = Alocar_vetor_real (ordemMatriz);
	dif = Alocar_vetor_real (ordemMatriz);
	vetor_filaAvaliacao = Alocar_vetor_real (ordemMatriz);
	rcv_disp = Alocar_vetor_int(noProcesses);
        rcv_cnt = Alocar_vetor_int(noProcesses);
	
	//Leitura da matriz e vetor
	Le_matriz_vetor(ordemMatriz,vetor,matriz,fp);

	//nao faz parte do algoritimo em si, apenas recolhemos valores iniciais pois está
	//na descricao do prloblema, é obrigatorio na saida do problema
	if(processId==0)
	{
		for(i=0;i<ordemMatriz;i++)
			vetor_filaAvaliacao[i]=matriz[filaAvaliacao-1][i];
		valor_filaAvaliacao=vetor[filaAvaliacao-1];
	}
	
	
	
	
	//define os vetores utilizados em MPI_Allgetherv
        if (resto>0) {
                rcv_disp[0] = 0;
                if (processId == 0) rcv_cnt[0] = nl;
                else rcv_cnt[0] = nl+resto;
                for ( i = 1; i < noProcesses; i++ ) {

                        rcv_disp[i] = (int)BLOCK_LOW(i,noProcesses,ordemMatriz-resto)+resto;
                        if (processId == 0) rcv_cnt[i] = nl-resto;
                        else rcv_cnt[i] = nl;

                }
        }
        else {

                for ( i = 0; i < noProcesses; i++ ) {

                        rcv_disp[i] = (int)BLOCK_LOW(i,noProcesses,ordemMatriz);
                        rcv_cnt[i] = nl;

                }

        }


	//divide pela diagonal, cada nó divide o seu correspondente
	k=0;
	for(i=posicao;i<(posicao+nl);i++)
	{
		aux=matriz[i][i];
		for(j=0;j<ordemMatriz;j++)
			matriz[i][j]=matriz[i][j]/aux;
		x[k]=vetor[i]/aux;
		k++;
	}

	

	//enche o buffer para a proxima interaco
        MPI_Allgatherv(x,nl,MPI_FLOAT,xg,rcv_cnt,rcv_disp, MPI_FLOAT,MPI_COMM_WORLD);
	
	
	omp_set_num_threads(NTHREAD);
	#pragma omp parallel for shared(xg,ordemMatriz) private(i) firstprivate(vetor)
	for(i=0;i<ordemMatriz;i++)
		vetor[i]=xg[i];

	
	//Demais Interações
	k=1;
	
	while(k<interacoes)
	{
		
		#pragma omp parallel for shared(xg,ordemMatriz) private(i) firstprivate(result_antes)
		for(i=0;i<ordemMatriz;i++)
			result_antes[i]=xg[i];

		
		#pragma omp parallel for shared(posicao,nl,ordemMatriz,matriz,vetor) private(i,j) firstprivate(aux,x)
		for(i=posicao;i<(posicao+nl);i++)
		{
			aux=0;
			for(j=0;j<ordemMatriz;j++)
				if(i != j)
					aux=aux+matriz[i][j]*result_antes[j];
			//Xdepois(i)=Vetor(i) - Soma(j=0 -> n)de{matriz[i,j]*Xantes(j)} para i<>j
			x[i-posicao]=vetor[i]-aux;
		}

		MPI_Allgatherv(x,nl,MPI_FLOAT,xg,rcv_cnt,rcv_disp, MPI_FLOAT,MPI_COMM_WORLD);

		
		//critério de parada por erro
		#pragma omp parallel for shared(ordemMatriz,xg,result_antes) private(i) firstprivate(dif)
		for(i=0;i<ordemMatriz;i++)
		{
			dif[i]=xg[i]-result_antes[i];
			if (dif[i] <= 0) dif[i] = (dif[i] * -1); //abs(dif[i])
		}

		
		mr=Numero_maximo(ordemMatriz,dif)/Numero_maximo(ordemMatriz,xg);
		if (mr <= 0) mr = (mr * -1); //abs(mr)
		if(mr<=erro)
			break;
		k++;
	}
	
	//impressao resultados
	if(processId==0)
	{

		aux=0;
		for(i=0;i<ordemMatriz;i++)
			aux=aux+vetor_filaAvaliacao[i]*xg[i];

		printf("Iterations: %d\n",k);
		printf("RowTest: %d => [%f] =? %f\n\n",filaAvaliacao,aux,valor_filaAvaliacao);
	}
	
	MPI_Finalize( );
	// libera memoria
	fclose(fp);
        fflush(stdout);
	matriz = Liberar_matriz_real (ordemMatriz, ordemMatriz, matriz);
	vetor = Liberar_vetor_real (ordemMatriz, vetor);
	result_pos = Liberar_vetor_real (ordemMatriz, result_pos);
	result_antes = Liberar_vetor_real (ordemMatriz, result_antes);
	dif = Liberar_vetor_real (ordemMatriz, dif);
	vetor_filaAvaliacao = Liberar_vetor_real (ordemMatriz, vetor_filaAvaliacao);
	

	return 0;
}

float **Alocar_matriz_real (int m, int n)
{
	float	**v;  // ponteiro para a matriz 
	int   	i;    // variavel auxiliar    

 	if (m < 1 || n < 1) // verifica parametros recebidos
	{ 
	     printf ("** Erro: Parametro invalido **\n");
	     return (NULL);
     	}

  	// aloca as linhas da matriz 
  	v = (float **) calloc (m, sizeof(float *));
  	if (v == NULL) 
	{
	     printf ("** Erro: Memoria Insuficiente **");
	     return (NULL);
     	}
	
  	// aloca as colunas da matriz
  	for ( i = 0; i < m; i++ ) 
	{
      		v[i] = (float*) calloc (n, sizeof(float));
     		if (v[i] == NULL) 
		{
         		printf ("** Erro: Memoria Insuficiente **");
         		return (NULL);
         	}
      	}
  	return (v); // retorna o ponteiro para a matriz
}

float **Liberar_matriz_real (int m, int n, float **v)
{
	int	i;  // variavel auxiliar

  	if (v == NULL) return (NULL);
  	if (m < 1 || n < 1) // verifica parametros recebidos
	{  
    		printf ("** Erro: Parametro invalido **\n");
     		return (v);
     	}
  	for (i=0; i<m; i++) free (v[i]);// libera as linhas da matriz 
  	free (v);      			// libera a matriz 
  	return (NULL);			// retorna um ponteiro nulo 
}

float *Alocar_vetor_real (int n)
{
	float	*v; // ponteiro para o vetor

	if (n < 1) // verifica parametros recebidos
	{  
	     printf ("** Erro: Parametro invalido **\n");
	     return (NULL);
     	}

  	// aloca o vetor
  	v = (float *) calloc (n+1, sizeof(float));
  	if (v == NULL) 
	{
	     printf ("** Erro: Memoria Insuficiente **");
	     return (NULL);
     	}
  	return (v);    /* retorna o ponteiro para o vetor */
}

int *Alocar_vetor_int (int n)
{
	int	*v; // ponteiro para o vetor

	if (n < 1) // verifica parametros recebidos
	{  
	     printf ("** Erro: Parametro invalido **\n");
	     return (NULL);
     	}

  	// aloca o vetor
  	v = (int *) calloc (n+1, sizeof(int));
  	if (v == NULL) 
	{
	     printf ("** Erro: Memoria Insuficiente **");
	     return (NULL);
     	}
  	return (v);    /* retorna o ponteiro para o vetor */
}

float *Liberar_vetor_real (int n, float *v)
{
	if (v == NULL) return (NULL);
  	if (n < 1) // verifica parametros recebidos
	{ 
     		printf ("** Erro: Parametro invalido **\n");
     		return (NULL);
     	}
  	free(v);        // libera o vetor
  	return (NULL);  // retorna o ponteiro
}

void Le_matriz_vetor(int n,float *a, float **b, FILE *fp)
{
	int	i,j; //manipulacao matriz e vetor

	for(i=0; i< n;i++)
		for(j=0;j<n;j++)
			fscanf(fp,"%f",&b[i][j]);

	for(i=0;i<n;i++)
		fscanf(fp,"%f",&a[i]);

}
float Numero_maximo(int n,float *v)
{
	int i;
	float saida;
	saida=v[0];
	for(i=1;i<n;i++)
		if(saida<v[i])
			saida=v[i];
	return saida;			
}
/*for(i=0; i< ordemMatriz;i++)
	{
		for(j=0;j<ordemMatriz;j++)
			printf("%f ",matriz[i][j]);
		printf("\n");
	}*/
